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Abstract. Moving puncture coordinates are commonly used in numerical simulations 
of black holes. Their properties for vacuum Schwarzschild black holes have been 
analyzed in a number of studies. The behavior of moving-puncture coordinates in 
spacetimes containing matter, however, is less well understood. In this paper we 
explore the behavior of these coordinates for Oppenheimer-Snyder collapse, i.e., the 
collapse of a uniform density, pressureless sphere of dust initially at rest to a black 
hole. Oppenheimer-Snyder collapse provides a stringent test of the singularity-avoiding 
properties of moving-puncture coordinates, since the singularity can form more quickly 
than it would for matter with pressure. Our results include analytical expressions for 
the matter density, lapse function, and mean curvature at early times, as well as 
interesting limits for later times. We also carry out numerical simulations to obtain 
the full solution and these show that, even in the absence of pressure, moving-puncture 
coordinates are able to avoid the singularity. At late times the geometry settles down 
to a trumpet slice of a vacuum black hole. 



1. Introduction 

Numerical simulations of black holes have recently experienced a dramatic breakthrough 
(see [H El |3] as well as numerous later publications). Many of these simulations now 
adopt some variation of the BSSN formulation jU [3] together with moving-puncture 
coordinates (2j |3] to handle the black hole singularities (see also [6] for a review). 

Moving-puncture coordinates adopt the 1+log slicing condition for the lapse 
function [7J and a gamma-driver condition for the shift vector [8]. The properties of 
moving-puncture coordinates have been studied in detail for Schwarzschild spacetimes 
[HI EH [HI H21 H3]- In particular, these studies show that dynamical evolutions of 
Schwarzschild black holes with moving-puncture coordinates settle down to a slice that 
ends on a sphere of finite areal radius. When displayed in an embedding diagram 
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(see, e.g., Figure 2 in [12]), these slices resemble a trumpet, which motivates the name 
"trumpet geometry" . The properties of these slices, in particular the fact that they do 
not extend to the singularity but instead end at finite areal radius, explain why these 
slices are well suited for dynamical simulations. 

Studying interactions of matter and black holes, which is important for many 
problems of astrophysical interest, requires relaxing the assumption of pure vacuum 
and including matter sources in the numerical simulations. While numerous simulations 
have already adopted moving-puncture coordinates in evolution calculations involving 
matter (see, e.g., [HI [151 HS1 E3 HH] ) , the properties of these coordinates in non- vacuum 
spacetimes are less well understood than in vacuum spacetimes. 

In a recent study, Thierfelder et al [19] analyzed the collapse of a spherically 
symmetric, unstable fluid star using two variations of moving-puncture coordinates (see 
their equation (2)). In the version that is similar to the gamma-driver condition most 
commonly used in binary black hole simulations (/is = 1 in their equation (2); see also 
|20j), the grid-points are pulled out of the star at late times during the collapse (see also 
[21] as well as the discussion at the end of Chapter 4 of [6]). All matter is effectively 
removed from the computational grid, and the simulation settles down to a trumpet 
geometry that becomes indistinguishable from that of a vacuum black hole. During this 
process the stellar interior becomes poorly resolved; this can make the code crash unless 
certain hydrodynamical values are prevented from taking unphysical values (see [19] for 
details). In an alternative version of the gamma-driver shift condition = a 2 in their 
equation (2)), the grid points are not pulled out of the star. Instead, the grid-points are 
"frozen" at late times and remain in the stellar interior; this approach leaves the matter 
well resolved even at late times. 

A classical, analytical solution describing the collapse of matter into a black hole is 
Oppenheimer-Snyder (OS) collapse of a homogeneous dust sphere [22]. This solution has 
served as a testbed for numerous numerical codes, and it has proven useful to transform 
this solution, originally derived in Gaussian normal coordinates, into coordinate systems 
that are more commonly used in numerical simulations (including maximal slicing [23J , 
polar slicing [21] and "observer-time" coordinates (22]). It is therefore of interest to 
study OS collapse in moving-puncture coordinates as well, which we do below. 

Moving-puncture coordinates form a rather complicated set of equations. It appears 
impossible to obtain a complete analytical solution to the OS collapse problem in these 
coordinates. However, we can obtain some partial analytical results and we do so here to 
illuminate our numerical integrations. Specifically, we find that at early times the matter 
density, lapse function, and mean curvature remain constant in space while varying in 
time. We also show that these functions remain spatially constant until a gauge mode 
has propagated from the stellar surface. Finally, we derive a lower limit for the lapse 
function at the center of the dust cloud. 

OS collapse poses another interesting question: can the lapse in moving-puncture 
coordinates collapse sufficiently fast for the slicing to avoid the black hole singularity, 
and to exhibit a transition to a trumpet geometry? As demonstrated by [19j, this is 
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Figure 1. Spacetime diagrams for OS collapse, for an initial areal radius of R = 5M. 
The dotted lines show trajectories of dust particles, which correspond to constant 
values of the radial coordinate \. The solid line represents a hypothetical surface of 
constant coordinate time t. 



the case for collapse of a fluid star with pressure. In the absence of pressure, however, 
as in OS collapse, the matter collapses more rapidly, leaving the lapse less time to drop 
to zero and avoid the singularity. In this sense, OS collapse poses a more stringent test 
than the collapse of a fluid star with pressure. We demonstrate numerically that, even 
for OS collapse, moving-puncture coordinates do avoid the singularity and allow for a 
transition to a trumpet geometry. 

This paper is organized as follows. In Section [2] we review the OS solution and 
introduce a transformation to general coordinates. In Section [3] we specialize to the 
1+log slicing condition. We show that the density, lapse, and mean curvature remain 
spatially constant at each point of the stellar interior until that point is reached by 
a gauge mode propagating from the stellar surface. In Section [4] we present the full 
numerical solution. We use this solution to probe the late-time behavior of OS collapse, 
demonstrating the transition to a trumpet geometry. We briefly summarize our results 
in Section |5j Throughout this paper we use geometrized units in which G = c = 1. 

2. Oppenheimer-Snyder collapse 

2.1. Solution in Gaussian normal coordinates 

OS collapse [22] describes the collapse from rest of a constant-density dust sphere to a 
black hole (see Figure [I] for a spacetime diagram, see also Section 1.4 in [U] for a review). 
The solution can be given analytically by matching a closed-Friedmann solution for 
the stellar interior to a Schwarzschild solution for the exterior. The interior solution, 
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expressed in Gaussian normal coordinates, is given by 

ds 2 = -dr 2 + a 2 (r) (dx 2 + sin 2 X d^ 2 ) , (1) 

where < % < xo is a radial coordinate that is comoving with each dust particle. The 
scale factor a(r) can be expressed parametrically as a function of the proper time as 
measured by each dust particle with the help of the conformal time rj, 

1 , 

a = -a m (l + cost]), 

\ (2) 
T = -a m (r] + smr]). 

Matching this solution to an exterior Schwarzschild solution shows that the initial scale 
factor a m and the maximum value xo of the radial coordinate \ are given by 



R 3 \ 1/2 



2M 



(3) 



and 



*» =sin WflbJ- (4) 

where Rq is the initial areal radius of the cloud, and M is the cloud's ADM mass. 

The rest-mass density po, as observed by an observer comoving with the matter, 
remains homogeneous (i.e. spatially constant) on each slice of constant r, and is given 
by 

| 3 - (5) 



Po{r) 



Po(0) V a ( r )/ 

Here the initial density at r = 0, po(0), is related to the initial areal radius Rq and the 
mass M by 

4:71 

M = — po(0)iJg. (6) 

We note that, for dust, p is equal to the total mass-energy density as observed by a 
comoving observer, so that po = u a UbT ab , where u a is the dust particles's four-velocity 
and T ab the stress-energy tensor. 

All particles reach the singularity with infinite density Po{t) when a = 0. According 
to equation <^2h this corresponds to a conformal time of i] = ir, or a proper time 

— = f — V /2 7T. (7) 

M \2MJ v ; 

In Gaussian normal coordinates, the mean curvature K = j^Ky, where is the 

spatial metric and the extrinsic curvature, is given by 

3 da , . 

K GN = — —. 8 
a dr 

We also note that the "Gaussian normal" normal vector Uq N , which is orthogonal on 
slices of constant proper time r, is aligned with the dust particles' four- velocity u a , 

n a GN = u a . (9) 
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2.2. General coordinate transformation 



We first consider a general coordinate transformation to a new time coordinate t and 
a new radial coordinate r. We have included a hypothetical surface of constant t in 
Figure [TJ We assume that we can express the old coordinates r and x as functions of 
the new coordinates, 



r = r 



(t,r) and x = x(t,r); 



(10) 



we will also assume that t = for r = and r = at the origin x = 0- Inserting these 
functions into the line element ([I]) results in the new line element 

ds 2 = -(f 2 - a 2 x 2 )dt 2 + 2(a 2 xx' - fr')drdt 



+(aY 2 



,12 



)dr 2 + a 2 sin x<^ 2 



11 



where we have defined the partial derivatives 



X 



dr 
di 

dx 
dt 



r = 



X' 



dr 
dr 

dx 
dr 



(12) 



Comparing the line element (11) with the 3 + 1 form of the line element 

ds 2 = -a 2 dt 2 + ^(dx 1 + /3 i dt)(dx j + ftdt) 
we first identify the spatial metric 

j • ( 2 12 12 2-2 2-2 ■ 2 n\ 

jij = diag(a x — t , a sm x, a sm x sin V) 
and its determinant 

I 2 12 >2\ 4 • 4 ■ 2 a 

7 — ( a X -t Jo sin X sm 



(13) 
(14) 
(15) 



We can then read off the radial, and only non-vanishing, component of the shift vector 



Pr 



2 • I 

a XX 



TT 



or 



a 2 XX' 



TT 



a 2 x' 2 — t' 2 



Finally, the lapse function is given by 

2 

a - 



^12 ^ 



XT') 



i\2 



(16) 



(17) 



3. Early stage of the collapse: analytical results 

3.1. The 1+log condition 

We now impose the 1+log slicing condition 

{dt - ftdi)a = -2aK. 
Here the mean curvature K = j^Kij can be expressed as 



;is) 
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We will also use the time evolution equation for K, 

(d t - 0%)K = -Y j DiD ja + aK l:j K ij + 4vra(p + 3), (20) 

where Di denotes the covariant derivative associated with the spatial metric 7^, 
is the extrinsic curvature, p = n a nbT ab is the matter density as observed by a normal 
observer, and S = 'j ab 'j a c T bc is the trace of the spatial stress, again as observed by a 
normal observer. 



We note that the slicing condition (18) can also be written as 

-^-a = n a d a a = -IK, (21) 

where n a is the normal vector on each slice of constant coordinate time t (which, in 
general, is different from Bq N ), and where r n is the proper time as observed by such 
a normal observer. This shows that 1+log slicing is invariant under spatial coordinate 
transformations, i.e. independent of the choice of the shift (see also [26]). In order to find 
the lapse a, we may therefore choose whatever shift is convenient, and we will consider 
two different choices below. 

We also observe that the pair of equations formed by the 1+log slicing condition 
(18) and the evolution equation for the mean curvature (20) can be combined to form 
a wave equation for the lapse. For simplicity we may assume vanishing shift. Taking a 
time derivative of equation (18) and using the result to eliminate the term d^K in (20) 
we obtain an inhomogeneous wave equation of the form 

- d 2 t a + 2«7 ij DiDja = (22) 

where the right hand side contains lower-order terms that do not affect the wave operator 
on the left-hand side. 



3.2. Spatially constant lapse 

We first show that if a, K and po are spatially constant in a region of a spatial slice, 
then they will remain spatially constant in the domain of dependence of this region. 
These quantities may then change in time, but not in space. This result can be verified 
directly from the equations; it is easier, however, to see this from the properties of OS 
collapse. 

According to equation ([5]), slices of constant po in OS collapse correspond to slices of 
constant proper time r (which is measured along the trajectories of dust particles). The 
region that we are considering therefore corresponds to a region of a slice of constant 
t. Moreover, the normal on these slices is aligned with the four-velocity of the dust 
particles, n a = u a , which implies p = po- According to equation (J8|, K is also constant 
on such a slice. We now assume a more general, but still spatially constant, lapse. 
Integrating forward to a new time slice of constant coordinate time t + dt, the proper 
time as measured along the slice's normal (and hence along the dust particle's trajectory) 
will increase by dr = adt. In a region where a is spatially constant, the new slice will 
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again correspond to a slice of constant r. Moreover, with K constant in this region, 



equation (21) shows that the lapse will again be constant on the new slice. 

In particular, this argument applies to early times, assuming the integration starts 
with a constant initial lapse ocq. On the other hand, this argument does not apply at 
the surface of the star. At the star's surface, the density p is discontinuous, and since 



this quantity appears on the right hand side of (20), the mean curvature will not remain 



spatially constant. We therefore expect that the surface will trigger a deviation that 



propagates into the star, satisfying the wave equation (22). The quantities p, K and 
a will remain spatially constant at any location in the star only until this deviation 
has arrived; i.e. in the domain of dependence of the initial stellar interior. We show a 
numerical demonstration of this behavior in Figure [2} 

Assuming a spatially constant lapse we can obtain its analytic value by choosing 



(3 r = 0. Equation (18) then simplifies and can be integrated immediately to yield 



a = a + ln(7/7o), (23) 



where the subscript denotes initial values at r = 0. Inserting the determinant (15) we 
find 

/ f a 2 y /2 - r /2 k 4 \ 

It is difficult to make further progress in general. If we assume, however, that the lapse 
has remained spatially constant, as it does at early times, then the above expression 
simplifies. While the lapse remains constant, slices of constant time t will coincide with 
slices of constant r, meaning that r' = and that n a = n^. For zero shift, coordinate 
observers are comoving with normal observers; this means that the new coordinates r 
are comoving with the fluid, which, in turn, are described by fixed values of the original 
radial coordinate %. For zero shift, and as long as the lapse remains constant, we 



therefore have x' — X01 an d the lapse (24) reduces to 

a = a + 61n(a/a m ) (as long as a remains spatially constant). (25) 
Alternatively, this expression can be derived by observing that, on slices of constant 



r, the mean curvature is given by (pft, so that for vanishing shift equation (18) can be 



integrated directly to yield (25). 



3.3. The lapse at the center 



We will now show that the expression on the right-hand side of Eq. (25) is a lower limit 
for the lapse function. To begin, we choose the shift vector in such a way that x = r at 
all times. In this case we have x — and x' = 1, so that the lapse reduces to 

a 2 f 2 

a 2 = (26) 

a 1 — r l 



and the shift to 
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Figure 2. Values for the lapse a for OS collapse with an initial areal radius of 
Ro = 5M. The left panel shows the values of the lapse at the center, a c , together 



with the lower limit «ll (see equation ( 30 1) . Note that the lapse (red) agrees with the 
analytical lower limit (green) at early times, but suddenly starts deviating from «ll at 
3.5M. The right panel shows profiles of the lapse as a function of areal 



a time r. 



gauge 



radius R at different times t, parametrized by the proper time r c at the center. Red, 
green, blue and magenta lines correspond to t/M = 0, 1.47, 2.38 and 3.54 respectively. 
At early times, the lapse remains spatially constant in a region decreasing in size 
around the center. At the time T gaugo this region disappears, and the lapse is no longer 



constant. 



Regularity at the origin demands /3 r = there, which implies r' = at r = (as long as 
f ^ 0). This means that slices of constant coordinate time t must be tangent to slices of 
constant r at the origin, as indicated in Figure [1} Inserting the lapse and shift, together 
with (19), into the 1+log slicing condition (18), and evaluating the resulting expression 
at the origin r = 0, we obtain 

1 da t"t , , 

f-6-— -2— =0. 28 
a dt a z 

This is a complicated equation for r, and it is again difficult to make progress, even 
after having restricted the analysis to the origin. However, we can derive a useful limit 
from the above expression. 

Let us assume that the lapse has never taken a local maximum at the center; in 
other words, we assume that the lapse has always been either spatially constant in a 
neighborhood of the center, or it has taken a local minimum. This condition holds at 
early times, while the lapse remains spatially constant, but even at later times this is 
a reasonable assumption for a slicing condition with "singularity-avoiding" properties; 
this is also what we find in our numerical simulations of OS collapse in moving-puncture 
coordinates. Under these conditions, coordinate slices of constant t will be "held back" 
at the center, as suggested in Figure [TJ and we will have 

t" > (29) 



at the center. We can then integrate equation (28) to obtain the limit for f > 
f(0) + 61n(a/a m ). Moreover, equation (26) demonstrates that a = f at the center; 
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hence 

ot c > oiq + 6 ln(a/a m ) = (as long as r" remains positive) (30) 

Not surprisingly, the lower limit a LL agrees with the value that we obtained under the 



assumption of spatially constant lapse in equation (25). 

In Figure [2] we show numerical results that confirm the above findings. These 
simulations were performed with a code that we describe in more detail in Section [4] 
below. In the left panel of Figure [2] we show both a c and «ll as a function of proper 
time. At early times, while the lapse remains spatially constant so that r" = 0, the 
two agree as expected. In the right panel of Figure [2] we show profiles of the lapse at 
different times t. These profiles demonstrate that the lapse indeed remains spatially 
constant in a region around the center. At a time that, for reasons that will become 
clear below, we will call the "gauge time" Tg auge , the central lapse suddenly starts to 
deviate from its lower limit. For Rq = 5M as in Figure [2] we have T gauge ~ 3.5M. From 
the right panel in Figure [2] we see that the region of spatially constant lapse disappears 
at the same time. 



3.4- The gauge time r. 



gauge 



We next analyze the gauge time T gauge . The numerical results in Figure [2] suggest that 
this time might correspond to the time at which the center has come into causal contact 
with the surface. 

We first compute the time at which a light ray, emitted from the surface of the 
star at t — 0, reaches the center. From equations ^ we have dr = adr), so that we 
can rewrite the Friedmann metric Q in terms of the conformal time i] in the conformal 
form 

ds 2 = a 2 (r) (-drj 2 + d^ + sin 2 X dtt 2 ) (31) 

(this motivates the name "conformal time" for 77). Evidently, an ingoing, radial light 
ray satisfies drj = —d\- A light ray sent out at conformal time rj = from the surface 
of the star at \ = Xo reaches the center at a conformal time rj = xo, or at a proper time 

^caus = ^a m (xo + sinxo)- (32) 
Inserting equations ^ and ^ we can write this in terms of i?o as 

Teaus / #0 \ 3/ % ( . -1 (2M\ 1 ' 2 /2M^ ^ 



[Mi) +U) ■ < 33 > 



M 

which, in the weak-field limit 2M/i? - > 0, yields the expected result r caus — >■ R . We 
have included this null characteristic as a dashed line in Figure |3j 

Evaluating (33) for R = 5M we find r caus = 5.21M (see Figure [3]), which is later 
than T gauge . This is not necessarily surprising, since the speed of light will only affect 
physical quantities. The lapse, however, is a gauge quantity, which may propagate 
at speeds faster than the speed of light. An analysis of the characteristic speeds for 
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Figure 3. A spacetime diagram of OS collapse for Rq = 5M as in Figure [T] except 
that here we include two characteristic curves originating from the stellar surface at 
t = 0. The dashed line marks a null characteristic, i.e. the trajectory of a light ray, 
while the solid line marks a superluminal gauge mode propagating into the star. As 
explained in the text, the density p, the lapse a and the mean curvature K remain 
spatially constant in the shaded region below the solid line, which has not yet come into 
"gauge" contact with the surface. Also included are four critical times at the center of 
the star, namely the time T gaugc at which the center comes into "gauge" contact with 
the surface, r caus when it comes into causal contact with the surface, Tij m j t when the 
1+log lapse drops to zero and halts any further evolution, and the proper time r s ; ng 
at which the singularity forms. 



the BSSN system in moving-puncture coordinates shows that some of them are indeed 
faster than the speed of light [2H EE]- For sufficiently weak gravitational fields the 
fastest characteristic, labelled %4 m [2H], has a proper speed, as measured by a normal 
observer whose four-velocity coincides with the normal on the slicing n a , of 



^char = -y2/at (34) 
(see Table I in |28j). 

At least heuristically, the existence of this characteristic speed can be motivated 
quite easily from the wave equation (22). From the principal part on the left hand side 
of this equation we find that radial modes propagate at coordinate speeds 

v™Z d = ±(2a 7 rr ) 1/2 - (35) 

Following the procedure described in [2H] these coordinate speeds can be converted into 
the proper speed (34). 

A normal observer comoving with the matter in the interior solution of OS collapse 
measures proper distances ad\ and proper times dr = adrj; such an observer therefore 
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observes this gauge mode propagating at a speed 



adx 
dr 



dx 
dq 



(36) 



We now assume that the lapse remains spatially constant in a region that has not yet 
come into "gauge" contact with the stellar surface; in this case the lapse is given by 
equation (25) and we can rewrite (36) as 



dx 



3 In 



1 + COS 7] 



-1/2 



drj, 



(37) 



where we have assumed ao — 1. It does not seem possible to evaluate this integral 
analytically. A simple numerical integration from xo to x — shows that, for R Q = 5M, 
this gauge mode arrives at the center at a conformal time ^ gaU ge = 0.456, or a proper 
time T gauge = 3.54M. This agrees well with the critical time identified earlier from our 
numerical simulations. We have included this characteristic as the solid line in Figure [3j 



4. Late stage of the collapse: numerical results 

4-1- Numerical Methods 

We performed numerical simulations of OS collapse in moving-puncture coordinates 
with two separate and independent codes, producing consistent results. 

One code was the 3D code developed by the Illinois group (see, e.g. [2H | |3"U | I3"T | [T7] ) . 
The code adopts the BSSN formalism (HE] to evolve the gravitational fields and a high- 
resolution shock-capturing scheme to evolve the relativistic hydrodynamics. In order to 
achieve sufficient numerical resolution the code adopts adaptive mesh refinement. For 
the results shown in this paper we used 8 levels of refinement, resulting in a resolution 
of SX/M = 0.0078125 on the finest level with the outer boundaries imposed at ±64M. 

Our other code is a ID code, originally developed for vacuum spacetimes in Ref. [28] . 
that uses the Cartoon method to impose spherical symmetry. We have supplemented 
this code with a simple finite-difference implementation of the equations of relativistic 
hydrodynamics assuming dust. All of the results shown in this Section were obtained 
with the 3D code, and have been confirmed with the ID code. 

We also show results for two different shift conditions. The first condition is the 
gamma-driver condition commonly used in moving-puncture simulations 

= f3 k d k f3 l + {d t T l ) rYls -f3 k d k T l -r,B l 

Here the auxiliary quantity B l has been introduced in order to turn an otherwise second- 
order in time equation into a first-order equation, the conformal connection functions 
f * = 7 /m f \ m are defined as contractions of the Christoffel symbols associated with the 
conformally related spatial metric 7^, and (d t T l ) T hs denotes their time derivative (see, 
e.g., Box 11.1 in [6]). Finally, we will adopt 77 = 2/M for all results shown here. 
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Figure 4. Conformal factor ip as a function of isotropic radius r for the initial dust 
sphere. The curves correspond to R — 5M, R = 3M and R = 2M; also included 
is the conformal factor i/j = 1 + M/(2r) for "wormhole" Schwarzschild data, which 
diverges for r — > 0. The dust curves match the wormhole curve at r = 3.94M (for 
i? = 5M), r = 1.87M (for R = 3M) and r = 0.5M (for i? = 2M). 



We also consider the alternative condition 
d t p = /i 5 f 1 - V (3 l + 



(39) 

3/4 [20] . Thierfelder e* 
a/ [19] considered two cases in their analysis of fluid stars: /is = 1, and fis = a 2 . We 



which is equivalent to the traditional condition (38) when 



also consider the condition (39) below; whenever we do, we adopt = a 2 . 

Initial data for our codes, which adopt cartesian coordinates, can be constructed by 
bringing the metric into isotropic coordinates. The exterior Schwarzschild metric then 
takes the familiar form 



dl 2 



1 + 



MV 



2rJ 



(dr 2 + r 2 dn 2 ), (40) 

where dl 2 denotes a spatial line element, and where r is an isotropic radius. In the 
interior, we can transform the Friedmann metric to isotropic coordinates by letting 
X = 2 arctan(cr) with the constant c chosen such that the conformal factors in the 
interior and exterior match at the star's surface. We then obtain the initial (spatial) 
line element 

dl 2 = ^(dr 2 + r 2 dn 2 ) 



with 




2M7.Ro) r R, 



2\ 

o 



1/2 



2rl + Mr 2 



r < r 



r > r 



(41) 



(42) 
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Figure 5. Profiles of the lapse a as a function of radius at different times t, as obtained 
with the shift condition (38). Red, green and yellow lines correspond to t/M = 0, 10, 
and 40, respectively. At late times, the lapse asymptotes to that the lapse of a trumpet 
slice of Schwarzschild in "stationary 1+log" slicing", which is included as a black line 
(see [H]). 



where ro = -Ro (l — M/Rq + \J\ — 2M/i?o) /2. The initial data also include Kij = 0. 
Profiles of the conformal factor if) for different values of Ro are shown in Figure |4j 

4.2. Results for R = 5M 

Results for the lapse a at early times have already been shown in Figure [2} All of 



these results have been obtained with the gauge condition (38), and demonstrate the 
existence of a region of decreasing radius around the origin in which these quantities 
remain spatially constant. We now turn to results for late times of the collapse. In this 
Section we will focus on R = 5M; we will discuss more general results below. 

In Figure [5] we show profiles of the lapse as a function of radius at different times t, as 



obtained with the shift condition (38). We also include the analytical lapse for trumpet 



slices of Schwarzschild in "stationary 1+log" slicing, which satisfies the equation 

2 2M C 2 e a 

a=1 ~if + ^- (43) 



with C = ( v / 2/16)(yi0 + 3) 3 / 2 e( 3 -^/ 2 M 2 w 1.24672M 2 (see, e.g., pj). Evidently, the 
lapse as obtained from OS collapse in moving-puncture coordinates asymptotes to that 
for a Schwarzschild black hole. This is not surprising for the stellar exterior, of course, 
but one might have expected that an interior region remains in which the density remains 
different from zero, and in which the lapse differs from that of a vacuum spacetime. As 
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Figure 6. Profiles of the density po as a function of areal radius R at different 



times t. The left panel shows result for the shift condition (381, while the right panel 



shows results for the shift condition ( 39 1 . For the former condition (left panel) , grid- 
stretching leads to all grid points being pushed out of the star at late times, leaving 
the star under-resolved, and making the evolution of the density very unreliable. To 
highlight this effect we marked individual gridpoints with crosses for the inner regions 
of our spacetime at t/M = 24. For the latter condition (right panel), grid points are 
frozen inside the star, so that it remains resolved even at late times. Red, green, blue 
and magenta lines correspond to t/M = 0, 14, 24 and 34 respectively. Note that there 
is no line corresponding to t.M — 34 in the left panel, as the matter is completely 
unresolved at this time. 



discussed in [T9] , however, the shift condition ( 38 ) pushes all grid points into the stellar 
exterior in a finite time, so that the stellar interior remains unresolved and removed 
from the grid. 

This effect is also visible in Figure [6j where we show profiles of the density po 
for different coordinate times t. In the left panel we again show results for the shift 



condition (38). At late times, the density drops to zero everywhere, indicating that all 



grid points having been pushed out of the stellar interior. In the right panel we show 



the same results for the shift condition (39). Here, grid points are frozen at late times, 



and remain in the stellar interior, which is well-resolved even at late times. 

In Figure [7] we show spacetime diagrams, similar to those in Fig. 6 of [12], to further 
illustrate the influence of the shift condition on the late-time trajectories of gridpoints. 
In each plot, the blue lines correspond to surfaces of constant areal radius, the green 
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Figure 7. Spacetime diagrams for OS collapse (compare Fig. 6 in [Ej). The left panel 



shows results for the shift condition (38), while the right panel shows results for the 



shift condition (39). Blue lines correspond to lines of constant areal radius. The green 



region marks the stellar interior. The red line denotes the apparent horizon which, as 
in Gaussian normal coordinates, first appears on the stellar surface. 



region marks the stellar interior, and the red line denotes the apparent horizon. In the 



left panel we show results for the shift condition (38), for which no gridpoints remain 



in the stellar interior at late times, while the right panel shows results for the shift 



condition (39), for which the stellar interior remains well resolved even at late times. 
The apparent horizon first appears when the stellar surface passes through R = 2M, as 
in Gaussian-normal coordinates (compare, e.g., Fig. 1.3 in [6]). 

In Figure M we show the central density function of time, as obtained with 



the shift condition (39). In the left panel we show the central density as a function of 



proper time, and compare with the analytical result ([5]). We see that the numerical 
evolution terminates at the proper time Ti imit (see also the spacetime diagram in Figure 
3]), when the lapse collapses to zero and prevents any further advance in proper time 
(recall dr = adt). In the right panel we show the central density as a function of 
coordinate time, which continues to advance even after the lapse has collapsed to zero. 
In this graph, the density levels off to the limiting value Po imit , corresponding to the 
value of the central density at proper time ri imit . 

Finally, we show profiles of the shift (3 l in Figure M As before, we show results 



obtained with the shift condition (38) in the left panel. Note that at late times, these 



lines terminate at a finite value of the areal radius R, corresponding to the limiting 
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Figure 8. The central density as a function of time, obtained with the shift 
condition ( 39 ) with ^5 = a 2 . The left panel shows the central density as a function 



of proper time r. Also included is the analytic result ([5]). The numerical solution 
terminates at the proper time Tii m j t , when the lapse collapses to zero, preventing any 
further advance in proper time dr — adt. The right panel shows the central density 
as a function of coordinate time t. In this graph the density levels off to the constant 
value p{) lmlt , indicating that, while coordinate time continues to advance, proper time 
does not. In each plot, the numerical value is shown in red, while in the left panel the 
analytic value is also shown in green. 



surface of the trumpet solution. In the right panel, showing results obtained with (39) 



the shift for small R drops to zero at late times, indicating that the corresponding 
gridpoints are "frozen" in the interior of the star. 



4.3. Effect of R 

In the previous Section we showed that for OS collapse with an initial areal radius of 
Ro = 5M, the 1+log slicing condition successfully halts the advance of the spatial slices 
before the singularity forms. In this Section we consider other values of Rq; in particular 
we evaluate whether or not, for smaller initial radii, the singularity forms too quickly 
for the lapse to halt the evolution. 

We have investigated this question numerically. Our simulations show that the 
1+log condition has no difficulty handling the collapse of arbitrarily small dust spheres, 
which, perhaps, is not surprising. Figure [4] shows the initial conformal factor ip as 




Figure 9. Profiles of the shift at different coordinate times t. The left panel shows 
values of the shift as obtained with the shift condition ( 38 1 , while the right panel shows 
results for the condition (39). Red, green, blue and magenta lines correspond to t/M = 
0, 14, 24 and 34, respectively. 



a function of isotropic radius r for different initial areal radii Rq/M, as well as the 
conformal factor for "wormhole" Schwarzschild initial data. As Rq decreases and 
approaches the limit 2M, the initial data approach that of a black hole. The black hole 
is the most stringent scenario, for which the lapse must collapse as quickly as possible 
to keep the slice from reaching the singularity. We also see this from the proper time ([7]) 
that elapses before the singularity is reached. For Rq = 2M we have r s ; ng = 7rM, which 
is the same time required for a normal observer at the center of a Kruskal diagram of 
a Schwarzschild black hole to reach the singularity. Numerous numerical simulations 
of the Schwarzschild spacetime, evolving "wormhole" initial data with moving-puncture 
coordinates, demonstrate that the lapse collapses sufficiently fast to avoid the singularity 
(e.g. [9j); apparently the same holds for OS collapse for any value of i?o- 

For R = 2M, the initial ball of dust is just inside its own Schwarzschild radius. In 
the spirit of black hole "stuffing" [32J |33l [31] , we can evolve the data with or without 
the matter fields. Omitting the matter terms in the black hole interior and evolving 
the initial data as if they were vacuum data leads to a violation of the constraints in 
the black hole interior; however, these constraint violations do not affect the black hole 
exterior, and disappear even in the interior after a short evolution of the data. 
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5. Summary 

We analyze OS collapse, describing the collapse of a homogeneous dust sphere initially 
at rest to a black hole, in moving-puncture coordinates. These coordinates have proven 
particularly robust for handling vacuum black holes in numerical simulations and this 
has motivated our interest in exploring how they behave in the presence of matter 
sources. 

We obtain the complete solution for OS collapse numerically and elucidate some 
of the key results analytically, especially for the early-time behavior. In particular, we 
show that if the initial lapse a® is chosen constant, then the lapse a, the density p and 
the mean curvature K will remain spatially constant in a region around the center of 
the star that decreases in size, and that disappears at a critical time i~ gauge . In this 
region, which is marked as the shaded region in Figure [3j the lapse is given by equation 



(25). The region is limited by a gauge mode that originates at the surface of the star 



and travels inwards at a speed that exceeds the speed of light (see equation (34)). We 
identify T gaugc with the time at which this characteristic, marked by the solid line in 
Figure [3j reaches the stellar center. At the center of the star, equation (30) serves as 
a lower limit on the lapse; the lapse agrees with this lower limit for times up to T gauge , 
but starts to deviate abruptly at that time. 



With the gamma-driver shift condition (38), grid points are pulled out of the stellar 



interior at late times, so that the numerical grid covers only the exterior of the star. The 
numerical solution asymptotes to a trumpet slice of a Schwarzschild spacetime. This 
confirms the findings of [19] and demonstrates that, even in the absence of pressure, 
the lapse collapses sufficiently fast to avoid the forming singularity, and to allow the 
transition to a trumpet geometry. We also confirm that with the alternative shift 



condition (39) with ps = en the coordinate grid points remain in the stellar interior and 



the star is well-resolved even at late times. 
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